home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / arith / ikarat.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  5.7 KB  |  238 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  ikarat.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16. #include <LEDA/impl/iint.h>
  17. #include <LEDA/impl/iloc.h>
  18.  
  19. #ifndef KARATSUBA_LIMIT2
  20. #define KARATSUBA_LIMIT2 15
  21. #endif
  22.  
  23. /* Formulae:
  24.     a == B^n a1 + a0;
  25.     b == B^n b1 + b0;
  26.     a b == (B^(2n) + B^n) a1 b1 + B^n (a1 - a0) (b0 - b1) + (B^n + 1) a0 b0;
  27.  
  28.     -----------------
  29. V    |     a1 b1     |
  30.     -----------------
  31.         -----------------
  32. IV        |     a1 b1    |
  33.         -----------------
  34.         -----------------
  35. III        ||a1-a0||b0-b1|    |
  36.         -----------------
  37.         -----------------
  38. II        |     a0 b0    |
  39.         -----------------
  40.             -----------------
  41. I            |     a0 b0    |
  42.             -----------------
  43.     ---------------------------------
  44.     |        prod        |
  45.     ---------------------------------
  46.  
  47.         |    |    |    |
  48.         n3    n2    n    0
  49.  
  50. */
  51.  
  52.  
  53. static  PLACE vecsub (/*diff, a, b, k*/);
  54. static  PLACE vecadd4carry (/*sum, a, b, c, d, carry, k*/);
  55. static  int vecadd3subcarry (/*accu, a, b, c, d, carry, k*/);
  56. static  PLACE vecaddPLACE (/*accu, a, d, k*/);
  57. static  PLACE vecaddscarry (/*accu, a, d, k*/);
  58.  
  59. void karatsuba_mult (prod, a, b, tmp, n2)
  60.     PLACE *prod, *a, *b, *tmp;
  61.     int n2;
  62.  
  63. /*     In:    a[n2], b[n2];    n2 == 2^k * n0, with n0 <= KARATSUBA_LIMIT2
  64.     Out:    prod[2*n2] = a[n2] * b[n2];
  65.     Local:    tmp[2*n2] for intermediate results
  66.     for n2>KARATSUBA_LIMIT2:    karatsuba-recursion with k-1
  67.     else:        standard multiplication
  68. */
  69. {    int sign, i, n, n3;
  70.     PLACE carry;
  71.     int subcarry;
  72.  
  73.     if (n2<=KARATSUBA_LIMIT2) {
  74.         /* use vecmuladd */
  75.         for ( i = 0; i < n2; i++)
  76.                 prod[i] = 0;
  77.             for ( i = 0; i < n2; i++) 
  78.                 prod[i+n2] = vecmuladd(&prod[i], a, b[i], n2);
  79.         return;
  80.     }
  81.  
  82.     n=n2/2;
  83.     n3=n2+n;
  84.  
  85.     karatsuba_mult(&prod[n2], a, b, prod, n);
  86.     /* prod:    |    a0*b0    |        |    */
  87.  
  88.     if (vecgt(a, &a[n], n)) {
  89.         sign=1;
  90.         vecsub(tmp, a, &a[n], n);    /* a0 - a1 */
  91.     } else {
  92.         sign=0;
  93.         vecsub(tmp, &a[n], a, n);    /* a1 - a0 */
  94.     }
  95.     /* tmp:        |    |    |    ||a0-a1||    */
  96.     if (vecgt(&b[n], b, n)) {
  97.         sign^=1;
  98.         vecsub(&tmp[n], &b[n], b, n);    /* b1 - b0 */
  99.     } else {
  100.         vecsub(&tmp[n], b, &b[n], n);    /* b0 - b1 */
  101.     }
  102.     /* tmp:        |    |    ||b1-b0|||a0-a1||    */
  103.  
  104.  
  105.     karatsuba_mult(&tmp[n2], tmp, &tmp[n], prod, n);
  106.     /* tmp:        ||b1-b0|*|a0-a1|||b1-b0|||a0-a1||    */
  107.  
  108.     karatsuba_mult(tmp, &a[n], &b[n], prod, n);
  109.     /* tmp:        ||b1-b0|*|a0-a1||    a1*b1    |    */
  110.  
  111.     /*
  112.             ---------------------------------
  113.     prod        |     a0 b0    |        |
  114.             ---------------------------------
  115.  
  116.             ---------------------------------
  117.     tmp        ||a1-a0||b0-b1|    |     a1 b1    |
  118.             ---------------------------------
  119.     */
  120.  
  121.     if (sign==0) {        /* Einfach addieren */
  122.         for (i=0; i<n; i++)
  123.         prod[i]=prod[n2+i];            /* I */
  124.         /* I + II + III + IV */
  125.         carry=vecadd4carry(&prod[n], 
  126.             &prod[n3], &prod[n2], &tmp[n2], tmp, 0, n);
  127.         /* carry + II + III + IV + V */
  128.         carry=vecadd4carry(&prod[n2], 
  129.             &prod[n3], &tmp[n3], &tmp[n], tmp, carry, n);
  130.         /* carry + V */
  131.         vecaddPLACE(&prod[n3], &tmp[n], carry, n);
  132.     } else {        /* III subtrahieren */
  133.         for (i=0; i<n; i++)
  134.         prod[i]=prod[n2+i];            /* I */
  135.         /* I + II + IV - III */
  136.         subcarry=vecadd3subcarry(&prod[n], 
  137.             &prod[n3], &prod[n2], tmp, &tmp[n2], 0, n);
  138.         /* carry + II + IV + V - III */
  139.         subcarry=vecadd3subcarry(&prod[n2], 
  140.             &prod[n3], &tmp[n], tmp, &tmp[n3], subcarry, n);
  141.         /* subcarry + V */
  142.         vecaddscarry (&prod[n3], &tmp[n], subcarry, n);
  143.     }
  144. }
  145.  
  146. static  PLACE vecsub (diff, a, b, k)
  147.         register PLACE *diff, *a, *b; 
  148.         register int k;
  149.         /* diff[k]=a[k]-b[k]; return carry (0 oder 1) */
  150. {       register PLACE carry=0;
  151.         for ( ; k>0; k--)
  152.         carry=PLACEsub(*a++, *b++, carry, diff++);
  153.         return (carry & 1);
  154. }
  155.  
  156. static  PLACE vecadd4carry (accu, a, b, c, d, carry, k)
  157.     PLACE *accu, *a, *b, *c, *d, carry;
  158.     int k;
  159.     /* accu[k]=a[k]+b[k]+c[k]+d[k]+carry; return carry; */
  160. {    PLACE ac, bc, cc, dc;
  161.     int i;
  162.     if (k==0)
  163.         return;
  164.     ac=PLACEadd(carry, *a, 0, accu);
  165.     bc=PLACEadd(*accu, *b, 0, accu);
  166.     cc=PLACEadd(*accu, *c, 0, accu);
  167.     dc=PLACEadd(*accu, *d, 0, accu);
  168.     for (i=1; i<k; i++) {
  169.         ac=PLACEadd(a[i], b[i], ac, &accu[i]);
  170.         bc=PLACEadd(accu[i], c[i], bc, &accu[i]);
  171.         cc=PLACEadd(accu[i], d[i], cc, &accu[i]);
  172.         dc=PLACEadd(accu[i], 0, dc, &accu[i]);
  173.     }
  174.     return ac+bc+cc+dc;
  175. }
  176.  
  177.  
  178. static  int vecadd3subcarry (accu, a, b, c, d, carry, k)
  179.     PLACE *accu, *a, *b, *c, *d;
  180.     int carry, k;
  181.     /* carry >= -1 */
  182.     /* accu[k]=a[k]+b[k]+c[k]-d[k]+carry; return carry; */
  183. {    PLACE ac=0, bc=0, cc=0, subc=0;
  184.     int i;
  185.     if (k==0)
  186.         return;
  187.     if (carry >= 0) {
  188.         ac = PLACEadd(carry, *a, 0, accu);
  189.         bc = PLACEadd(*accu, *b, 0, accu);
  190.         cc = PLACEadd(*accu, *c, 0, accu);
  191.         subc = PLACEsub(*accu, *d, 0, accu);
  192.     } else {
  193.         ac = 0;
  194.         bc = PLACEadd(*a, *b, 0, accu);
  195.         cc = PLACEadd(*accu, *c, 0, accu);
  196.         subc = PLACEsub(*accu, *d, 1, accu);
  197.     }
  198.     for (i=1; i<k; i++) {
  199.         ac=PLACEadd(a[i], b[i], ac, &accu[i]);
  200.         bc=PLACEadd(accu[i], c[i], bc, &accu[i]);
  201.         cc=PLACEadd(accu[i], 0, cc, &accu[i]);
  202.         subc=PLACEsub(accu[i], d[i], subc, &accu[i]);
  203.     }
  204.     return ac+bc+cc-subc;
  205. }
  206.  
  207. static  PLACE vecaddPLACE (accu, a, d, k)
  208.     PLACE *accu, *a, d;
  209.     int k;
  210.     /* accu[k]=a[k]+d; return carry; */
  211. {    PLACE carry=d;
  212.     for ( ; k>0; k-- )
  213.         carry=PLACEadd(*a++, carry, 0, accu++);
  214.     return carry;
  215. }
  216. /* ACHTUNG:    k==0  ==>  carry==d */
  217.  
  218. static  PLACE vecaddscarry (accu, a, scarry, k)
  219.     PLACE *accu, *a;
  220.     int scarry, k;
  221.     /* RADIX > scarry>=-1 */
  222.     /* accu[k]=a[k]+scarry; return carry; */
  223. {    PLACE carry;
  224.     if (scarry >= 0) {
  225.         carry = scarry;
  226.         for ( ; k>0; k-- )
  227.             carry=PLACEadd(*a++, carry, 0, accu++);
  228.         return carry;
  229.     } else {
  230.         carry = 1;
  231.         for ( ; k>0; k-- )
  232.             carry=PLACEsub(*a++, 0, carry, accu++);
  233.         return carry;
  234.     }
  235. }
  236. /* ACHTUNG:    k==0  ==>  carry==d */
  237.  
  238.